% function im = preconditionedAdjointPPFT2D(pp1,pp2) 
%
% Computes the preconditioned adjoint of the 2-D pseudo-polar Fourier transform.
%
% pp    2-D array of size (2n+1)x(n+1) containing the 2-D
%       pseudo-polar Fourier transform
%
% See also PPFT3, adjPPFT3, optimizedAdjPPFT3.
%
% XXX NOTE!!!!
% Compare to precondAdjPPFT3 in the folder Radon3. Note that n was
% replaced with n+1. This does not affect anything since n is even (since
% we take 3*n/2 and expect to get an integer), and so floor((n+1)/2) is
% exactly the same as floor(n/2). (July 2nd, 2007)
%
% Yoel Shkolnisky 28/02/03


function im = preconditionedAdjointPPFT2D(pp1,pp2,precision)


n=verify2DPPInput(pp1,pp2);
n=n(1);
m=2*n+1;
alpha = 2*(n+1)/(n*m);

%Compute the adjoint of PP1
tmp = zeros(2*n+1,n,precision);
for k=-n:n
   u = fliplr(pp1(toUnaliasedIdx(k,2*n+1),:));
   v = mult(k,alpha,n).*cfrft(u,-k*alpha,precision);
   tmp(toUnaliasedIdx(k,2*n+1),:) = v(1:n);
end

tmp = m*icfftd(tmp,1,precision);
adjpp1 = flipud(tmp(n/2+1:3*n/2,:));

%Compute the adjoint of PP2

tmp = zeros(2*n+1,n,precision);
for k=-n:n
   u = fliplr(pp2(toUnaliasedIdx(k,2*n+1),:));
   v = mult(k,alpha,n).*cfrft(u,-k*alpha,precision);
   tmp(toUnaliasedIdx(k,2*n+1),:) = v(1:n);
end
% To follow the code in adjPPFT we should have transposed each row before we assign it to tmp 
% (and creating an array of size nx(2n+1)). Then we had to apply cfft along the rows.
% To save operations, we assign the vector v to tmp without transpose, apply cfft along columns
% and then transpose the entire matrix at once.

tmp = m*icfftd(tmp,1,precision);
tmp = tmp.';
adjpp2 = flipud(tmp(:,n/2+1:3*n/2,:));

%Combine both adjoints
im = adjpp1+adjpp2;


function v = mult(k,alpha,n)
% Compute preconditioning factor for row k
if k==0
    v = 1/((2*n+1)^2);
else
    v = abs(k*alpha);
end